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Abstract. We analyze how the process of dust settling affects the spectral energy distribution and optical appear- 
ance of protoplanetary disks. Using simple analytic estimates on the one hand, and detailed l+l-D models on the 
other hand, we show that, while the time scale for settling down to the equator may exceed the life time of the 
disk, it takes much less time for even small grains of 0.1 /im to settle down to a few pressure scale heights. This is 
often well below the original location of the disk's photosphere, and the disk therefore becomes effectively 'flatter'. 
If turbulent stirring is included, a steady state solution can be found, which is typically reached after a few x 10^ 
years. In this state, the downward settling motion of the dust is balanced by vertical stirring. Dependent on the 
strength of the turbulence, the shape of the disk in such a steady state can be either fully flaring, or flaring only 
up to a certain radius and self-shadowed beyond that radius. These geometries are simila r to the geometries that 
were found for disks around Herbig Ae/Be stars in our previous papers fDuUemond l20o3 : Dullemond & Dominik 
A&A in press, henceforth DD04). In those papers, however, the reason for a disk to turn self-shadowed was by loss 
of optical depth through dust grain growth. Here we show that dust settling can achieve a similar effect without 
loss of vertical optical depth, although the self-shadowing in this case only affects the outer regions of the disk, 
while in DD04 the entire disk outside of the puffed-up inner rim was shadowed. In reality it is likely that both 
grain growth and grain settling act simultaneously. The spectral energy distributions of such self-shadowed — or 
partly self-shadowed — disks have a relatively weak far- infrared excess (in comparison to flaring disks). We show 
here that, when dust settling is the cause of self-shadowing, these self-shadowed regions of the disk are also very 
weak in resolved images of scattered light. A reduction in the brightness was already predicted in DD04, but we 
find that dust settling is far more efficient than grain growth at dimming the scattered light image of the disk. 
Settling is also efficient in steepening the spectral energy distribution at mid- to far-infrared wavelengths. From 
the calculations with compact grains it follows that, after about 10^ years, most disks should be self-shadowed. 
The fact that some older disks are still observed with the characteristics of flaring disks therefore seems somewhat 
inconsistent with the time scales predicted by the settling model based on compact grains. This suggests that 
perhaps even the small grains (< 0.1/xm) have a porous or fractal structure, slowing down the settling. Or it could 
mean that the different geometries of observed disks is merely a reflection of the turbulent state of these disks. 
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1. Introduction 

The dusty circumstellar disks surrounding most T Tauri 
stars and Herbig Ae/Be stars are believed to be the birth- 
place of planetary systems. These systems are therefore 
the object of intense study, both observationally and theo- 
retically. Much of what we know of these systems has been 
derived indirectly from their spectral energy distributions 
and infrared spectra. In the last few years the amount of 
observational information on spatial structure of these ob- 
jects has increased dramatically. The outer disk regions 



(typically outside of 50 AU) can be probed by direct 
imaging, both at optical/near- i nfrare d wavelengt hs wit h 



199£), 



19971) . 



HST/STIS (e.g. Augereau et al. l200H Grady et al. 
and in the sub-mm regime (Mannings & Sargent 
and with near- and mid-infrared interferometry t he in- 
ner disk structures can b e probed (e.g. Eisner et al. 120031: 
Millan-Gabet et al. 1200X1: Leinert et al. A&A subm.). 

Models of protoplanetary disks are increasingly suc- 
cessful at accounting for much of the observed proper- 
ties. For instance, they can explain why their spectral 
energy distributions (SEDs ) are generally rather flat in 
vF^ (Kenyon & Hartmann Il987j) . why dust featur es are 
almost always obser ved in emission (Calvet et al. Il99lt 
Chiang & Goldreich|l993), what the meaning is of the 



2 



DuUemond & Dominik: Dust settling in protoplanetary disks 



near- infrar ed bump in the SEDs of Herbig Ae st ars (Natta 
et al. l200lt DuUemond, Dominik & Natta l200l|) . and why 
some sources hav e stro ng and some have weak far-IR ex- 
cess (Duhemond '2002'; DuUemond & Dominik A&A in 
press, henceforth DD04). 

For the overaU disk properties, it is often sufficient to 
work with a disk model in which dust and gas are weU 
mixed and coupled. However, there are indications that 
dust growth, thermal processing, and dust-gas separation 
can play an important role in disks. The different overall 
shapes of the SEDs of Herbig AeBe stars (the "grou p I" 
and "group H" spectra identified by Meeus et al. l200l|) can 
be understood by a reduction in opacity in the o uter d isk 
via grain growth (DD04; see also D'Alessio et al. l200ll) or 
a direct reduction of the d isk height by settling of small 
grams (Chiang et al. l200lh . The high submm fluxes ob- 
served in disks with relatively low far-IR fluxes can be 
explained by a component of large grains in the disk in- 
terior which do not contribute to the opacity high above 
the midplane but contribute only as a co ld component 
to the submm emission (Natta et al. 120011: DD04). Both 
grain growth and grain settling is required to understand 
the presence of these grains. Moreover, the infrared spec- 
tra taken from the ground and from space have revealed 
structure in the infrared dust emission features that are 
indicative of dus t grain grow th and thermal proc essing 
(Bouwman l200l[ van Boekel l200.lt Hon da et al. 120031: 
Meeus et al. l2003| Przygodda et al. l2003(l . 

Also from a theoretical point of view, grain process- 
ing and dust-grain separation are expected to occur. In 
fact, a central role in our current understanding of planet 
formation is pla yed by th e forma.tion of a dust subdis k 
within the disk l)Safronovl Il969l: iDubruUe et"all Il995l) . 
Dust grains tend to settle towards the midplane, unless 
this process is countered by strong vertical stirring in the 
disk. The formation of such a subdisk with enhanced dust 
density greatly influences the timescales in which dust 
grains may grow to planetesimals and eventually plan- 
ets. For the detailed structure of the dust subdisk, pa- 
rameters like the remaining accretion rate, local turbu- 
lence. radial mixing processes, etc. play an d important role 
l|Dubrulle et all 1199,4 ICuzzi et alil200.ll) . Unfortunately, 
the processes that generate turbulence in protoplanetary 
disks are not well understood, and mixing/stirring is gen- 
erally treated in a simple parameterized manner. 

One way to approach this problem is to search for ob- 
servational signs of dust settling in disks. Observations of 
disks, at least in the optical and infrared, mostly probe the 
upper disk layers. If settling is active in disks, the upper 
layers of a disk will be influenced and may show observ- 
able changes. While there are many studies into the effects 
of settling onto the dust distribution relevant for planet 
formation (see e.g. Takeuchi & Lin 2_0_Q 2, 2003), there are 
curiously few studies which discuss the observational con- 
sequences of dust settling in a quantitative way. A st udy of 
this kind was presented by Miyake and Nakagawa |l99^ 
hereafter MN95). They studied the effects of grain set- 
tling and accretion luminosity on the SED of T Tauri 



disks in order to explain the infrared and IRAS colors of 
these stars. However, their study is limited in several ways. 
MN95 did not consider turbulent stirring at all. The basic 
assumption of their model was that turbulence has com- 
pletely ceased and dust settling can proceed undisturbed. 
The dust was assumed to be confined within one pressure 
scale height, and the time scales for changing the surface 
height in the disk were computed using the midplane den- 
sity. These time scales turn out to be between 10^ and 10^ 
years, and therefore comparable to the lifetime of proto- 
planetary disks. However, the timescales for settling are 
much shorter in the higher layers of the disk. Moreover, 
we will show that turbulence prevents the dust from set- 
tling down below a certain height, often much larger than 
the pressure scale height, and that turbulence will there- 
fore determine the final shape of the photosphcric surface 
of the disk. 

In this paper we investigate the effects of dust set- 
tling in the presence of turbulent stirring on the observ- 
able quantities of protoplanetary disks, in particular their 
SEDs and their images in scattered light. To do this we 
start with a well-mixed disk model for a T Tauri star and 
compute how the process of dust settling proceeds. In all 
stages of this process we compute observable quantities 
and see how they change. In doing this we hope to find 
typical observable signatures of dust settling which can be 
used to analyze protoplanetary disks. 

The paper is organized as follows: In section |21 we use 
simple analytical estimates to derive the timescales for 
settling and stirring. In section |31 we solve the equations 
for settling and stirring in a one-dimensional slab, a verti- 
cal cut through the disk and show how quickly the upper 
layers of the disk are depleted and to what depth. In sec- 
tion H we run a similar calculation at all locations in the 
disk and compute the changes settling causes to both the 
SED and images in scattered light. 

2. Settling, vertical stirring and accretion 

The equations for settling, vertical stirring and accretion 
have been extensively discussed in the literature. But, for 
the sake of clarity, and to define our notation, we include 
a brief discussion of the entire set of equations used for 
the current study. 

Dust grains in a protoplanetary disk experiences a 
gravitational force towards the star. If it were not for 
the gas in the disk, the grains would move on inclined 
orbits. In a co-rotating frame this means they would os- 
cillate about the midplane with the Keplerian frequency. 
However, small grains have a large surface-to-mass ratio, 
and therefore experience strong drag forces with the gas. 
Instead of oscillating about the midplane, these grains 
slowly settle toward the midplane on a time scale of a 
few hundred thousand years. This time scale depends on 
the surface density of the gas and on the surface-to-mass 
ratio of the grain. Large compact grains have a much 
smaller surface-to-mass ratio and therefore tend to set- 
tle on a much shorter time scale. Settling can, however, 
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only take place when turbulence in the disk does not stir 
up the grains above the midplane again. Usually this hap- 
pens up to a certain height Zgctt above the midplane while 
above this height the settling can proceed without prob- 
lems. This means that after a certain settling time, the 
disk above Zgctt is more or less devoid of dust grains of 
this particular size, while below it the vertical stirring has 
more or less equalized the abundance of dust grains of this 
size. 

In this section we derive the equations for these pro- 
cesses, and estimate the height Zgctt and the time scale 
that it takes for the dust to settle down to this height. 

2.1. Settling 

The interaction between a dust grain and the surrounding 
gas can be described by the friction time which is the time 
scale in which a dust grain responds to the motion of the 
gas. 



^fric 



mv 



(1) 



where m is the mass of the dust particle and v the rela- 
tive velocity between gas and dust. The friction force i^tric 
is a function of the drift velocity v, the grain properties 
and the gas density and temperature. In general, different 
regimes must be distinguished depending on the hydrody- 
namic regime appropriate for a given situation. However, 
for grains smaller than typically 1 cm, the Epstein regime 
is applic able throughout a protoplanetary disks (Cuzzi et 
so we restrict ourselves for the current paper to 
this regime. Then, under the generally valid assumption 
of subsonic drift velocities, the friction force is given by 



ric — -^PCrVCs 



(2) 



where p is the gas density, a is the collisional cross section 
of the dust grain and Cs is the isothermal sound velocity 



= (3) 

y /^gas"lp 

Tgas is the gas temperature, k the Boltzman constant, /Xgas 
the mean molecular weight in the gas and nip the proton 
mass. In the Epstein regime, the friction time can be writ- 
ten as 

3ml 

ifric = 7 . 4) 

4: a PCs 

The friction time depends upon the dust grain properties 
through the mass-to-crosssection ratio m/cr. For compact 
spherical grains with grain radius a and specific density 
Pd, raja = Ap^a/S, and the friction time turns into the fa- 
miliar form ifric = pdo,/ pCs- However, through much of this 
paper we will continue to write m/a since this expression 
is much more general and allows for a direct generalization 
towards 2-dimensional grains (like polycyclic aromatic hy- 
drocarbons (PAHs)) and porous or fractal grains. 

We now consider a circumstellar disk around a star 
with mass M^,. We assume the disk to be isothermal 



in the vertical direction, and we make the thin-disk- 
approximation: z ^ i?. At a distance R from the star, 
the vertical density distribution is then given by 



p{z) = 



2ttH, 



■ exp 



(5) 



where E is the surface density in the disk. Hp the pressure 
scale height and z the height above the disk. The pressure 
scale height is given by 



(6) 



where Q,k = \/GAh/ denotes the local Kepler fre- 
quency. A dust grain located at {R, z) is subject to a ver- 
tical gravitational force 



grav 



= — milT. 



(7) 



and, in the absence of turbulent stirring, settles towards 
the midplane at a settling speed 



Wsctt 



Apcs 



(8) 



which follows from 



grav fric ■ 

The settling speed decreases as one comes closer to the 
equator, because the gas density increases and the vertical 
component of the gravity decreases. The time for a grain 
to settle to a certain height z is therefore dominated by 
the local settling speed at height z 



tfsr.tt — 



Wsctt 



4 0- pCs 



(9) 



Using the vertical density distribution in Eq.©, we find 



3^/27^TO f^K 



exp 



2H^ 



(10) 



The time for a grain to settle towards the midplane 
of a disk is therefore approximately given by tgctt = 
(4/3-\/27r)(cr/m)(E/r2K). However, it is important to 
know, that in the higher layers of the disk, the timescales 
are very different. Typical flaring disk models intercept 
stellar radiation at a height z ~ 4iJp. The local settling 
time scales in this region differ from the time scales close 
to the midplane by a factor w 1/3000. Therefore, dust 
settling will deplete the upper layers of a disk surface much 
faster than the actual formation time of the dust subdisk 
close to the mid plane (see also Chiang & Goldreich ll997t 
Takeuchi & Lin l2002() . 

2.2. Vertical stirring 

Vertical stirring occurs when the disk is turbulent. The 
turbulent eddies of the gas transport the dust grains up 
and down in a random manner, thus constituting a diffu- 
sion process for the dust. A detailed formulation of the the- 
ory behind this tur bulent diffusion of dust grains is given 
by e.g. Volk et al. l)l980(l : Weidenschilling l|l997(l : Cuzzi, 
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Dobrovolskis & Champney ^1993); DubruUe, Morfill & 
Sterzik ( 199S) , and most recently by Schrapler & Henning 
(A&A subm.). 

The diffusion coefficient for gas molecules Dq is defined 
using the standard a prescription: 



Dq = aCsHp 



(11) 



But dust particles may not always couple perfectly to the 
gas. The actual diffusion coefficient for dust particles of a 
certain size is: 

Do 



D = 



(12) 



where Sc is the Schmidt number. For perfectly coupled 
particles (i.e. infinitely small grains) the Schmidt number 
is unity. The decoupling from the turbulence is described 
by a Schmidt number that is larger than unity. Various 
expressions for the Schmidt number have appeared in the 
literature. We adopt: 



Sc = (1 + St) , 



(13) 



where St is the Stokes number. The Stokes number St 
is defined as the ratio of the friction time tftic and the 
eddy turn-over time icdd- The latter depends not only 
on the turbulence parameter a, but also on the typical 
velocity of the largest eddies. This eddy velocity may de- 
pend on the type of turbulence (e.g. shear turbulence, con- 
vective turbulence, magneto-rotational turbulence), and it 
is not clearly agreed in the literature what it should be. 
Therefore we parameterize it in the following way: 



Vedd = oflCs , 



(14) 



where v^dd is the average velocity of the largest eddies and 
q a 'turbulence parameter' between and 1. Schrapler & 
Henning take q = 1/2 in their calculations. The typical 
largest eddy size Zedd can be found by the definition that 
Dq = ^cdd^cdd- One obtaius Zodd = a^^'iJp. This leads 
us to an expression for the eddy turn-over time iodd = 
ledd/vedd = Q!^~^''/f^K- The Stokcs numbcr then becomes: 



St = ^^"'^ = 3 m 17k^2<;-i 
todd 4 a pCs 



(15) 



In this paper we will take q — 1/2, following Schrapler & 
Henning. 

The turbulent a parameter used here is proportional to 
the ttaccr of the accretion process in the disk. Traditionally 
this is denoted as: 



Pt: 



(16) 



where Pt is the Prandtl numbcr, which is a dimcnsionless 
constant of order unity. In this paper we assume it to be 
Pt = 1 for the vertical turbulence. 



2.3. Equilibrium between stirring and settling 

High above the midplanc the densities of the gas are so 
low that turbulence is not able to effectively mix the grains 
upward. Instead the grains settle to the midplanc unhin- 
dered with a velocity given by Eq. ||SJ). At some height z 
above the equator the densities become high enough for 
the turbulent eddies to start to counteract the settling. 
This height can be estimated by comparing the time scale 
for settling igott to the time scale for diffusion tsta = /D 
which obeys 



is 



Sc 



(17) 



By equating igtir = ^^sott (where ^ is a factor described 
below) one obtains the following equation in z: 



Hi 



exp 



2m 



3\/27r m Sc(z) 
4 



(18) 



By taking ^ = 1 and solving for z (taking into account that 
Sc is a function of z), one obtains the height Zgctt/ R down 
to which the grains of a typical size and mass tend to set- 
tle. Below this height the turbulent diffusion manages to 
sweep up the grains and keep them afloat. The abundance 
for grains of this size is almost constant below this height. 
Above this height the abundance decreases because dust 
can settle down to lower altitude. It should be noted, how- 
ever, that the transition from fully mixed to fully depleted 
is not extremely sharp. As will be shown in Section |3| by 
solving Ea. (|18|l with ^ = 100 one obtains a good estimate 
of the 'depletion height' zjcpi above which the turbulence 
is so ineffective that dust grains of this size have depleted 
almost entirely. The transition from well-mixed to dust- 
free takes place between Zgott and Zdopi- Incidently, it is 
interesting to note that for grains smaller than 100 /im the 
Schmidt number Sc is almost always close to 1 at Zdcpi and 
below. The decoupling of the grains from the turbulence 
therefore has only marginal effect on the solutions. 

2.4. Steady state analytic estimates 

Based on the above equations one can now make figures 
of these quantities for a particular disk model. Consider 
a T Tauri star of mass Af* = 0.5AIq, radius i?* — 2.5Rq 
and effective temperature T* = 4000K. Surround it with 
a disk with a gas-|-dust surface density S(i?) given by 



S(-R) = So 



R 
1 AU 



(19) 



with p ~ —1.5 and Sq = 400g/cm^, with an inner radius 
of i?in = 0.1 AU and and out radius i?out = 400AU. The 
disk mass, for this configuration, is Afdisk = 0.011 Mq. We 
will assume, for the time being, that the disk is vertically 
isothermal, with a temperature T(R) that is determined 
by the irradiation of the disk by the central star. We ig- 
nore viscous heating, but we discuss its effect in Subsection 
15.21 In order to keep things simple and analytic, we as- 
sume that the flaring angle (p (see Chiang & Goldreich 
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I1997I for the definition) is a constant, and has a value of 
(f) = 0.05, which is not too far from realistic values. By 
equating the irradiation flux Fin. = (/)L*/(47ri?^) to the 
cooling flux Fcooi — o-sT{R)'^ one can solve directly for 
T{R), obtaining: 

T{R) = 4>^/^^T,. (20) 

We ignored geometric effects of order z/ R. Given the tem- 
perature T(R), the vertical density structure is then given 
by Eqs.(ini El ©■ Note that this disk structure is highly 
simplified, and in particular the assumption that </) is a 
positive constant is often untrue. But for the purpose of 
this paper it is sufficient. Throughout this paper the dust 
specific weight is taken to be 3.6 g/cm'^ and the mean 
molecular mass is taken /i = 2.3. 

In Fig. n the resulting Zdcpi/-R and the time scale at 
which the equilibrium state is reached are shown for the 
case of a = 0.01 and q — 1/2. In Fig.|5]the same is shown, 
but this time for a = 0.0001. Shown in the left panels 
in both figures is a grey line, which indicates the initial 
height of the r = 1 surface for radially outward moving 
stellar photons before the grains have started to settle. 
It is interesting to observe that for the a = 0.0001 case 
the height z^cpi/R below which the grain settle (even the 
0.1/im grains) is below this initial r = 1 surface height. 
Since dust grains carry virtually all of the opacity of the 
disk, this means that when the grains settle, the photo- 
spheric surface height Hg of the disk is reduced, i.e. the 
optical appearance of the disk becomes flatter. Also the 
infrared emission from the disk is reduced as the disk then 
captures less stellar flux. 

The steady state is always reached on a shorter time 
scale than the accretion time scale (taccr = R^ /v). It is 
shorter at small radii than at larger radii, but beyond 
a certain radius it stays constant. This is because these 
grains have settled below one pressure scale height, and 
for E IX. R~^-^ the time scale for a grain to settle down to 
one pressure scale height or l ower is independent of radius 
ijMivake fc Nakaeawal llQQSll . The settHng time scale for 
0.1 /im grains at a few hundred AU 10^ year) is not 
much smaller than the typical life time of protoplanetary 
disks (~ 10^ • • • 10'^ year). 

An important feature of the steady state solutions is 
that the Zdepi{R)/R does not always form a flaring shape. 
Beyond a certain radius iiturn(a) (different for different 
grain sizes), the value of d{zdcpi{R) / R) / dR becomes neg- 
ative. If the disk has only one grain size, and assuming that 
Zdcpi is a good estimate of the resulting surface height Hs 
of the disk, this would imply that the disk has an Hs/R 
that is not a monotonically increasing function of R. This 
is a strong indication that self- shadowing effects play a role 
in such disks: the outer disk has a lower Hs/R than the 
intermediate regions and therefore resides in the shadow 
of these intermediate disk regions (note that if one defines 
Hs to be the location of the r = 1 surface for radially- 
outward moving photons, a lower Hs/R at larger radii is 
meaningless because of the shadowing; see Fig. EJ. For 
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Fig. 3. The time scale with which grains settle below the 
T — 1 surface. The solid grey curve is the viscous (ac- 
cretion) time scale for a = 0.01. The dotted grey curve 
is the Kepler time scale tkcp — 1/^k- Since the settling 
cannot happen on a time scale smaller than this tkepi we 
have added ikcp to the plotted time scale to ensure this 
minimum. 



the case of a = 0.01 and 0.1 ^m grains this effect is small. 
One must go to 10 or 100 /.tm grains to see this effect. 
For a = 0.0001 this effect already takes place with 0.1 /.tm 
grains. 

In the calculations shown above the Stokes number 
remains mostly well below unity for grain sizes up to 100 
micron. However, if we take g = 0, the Stokes number 
becomes significantly larger than unity beyond the point 
where d{zdcpiiR) / R) / dR becomes negative. This increases 
the shadowing effect since the Zdopi(-R) drops more quickly 
as a function of R than for the case q — 1/2. The radius 
R turn of this turn-over point, however, does not change 
considerably. 

It should be noted that the time scale for the settling 
can be deceptive. It is the time scale it takes to reach 
the stationary state: the equilibrium between turbulence 
and settling. But the very upper parts of the disk get 
depleted from dust more quickly, since the dust settling 
velocity increases with the decreasing density. The time 
scale for settling of 0.1/im grains in the very outer regions 
of the disk is therefore an upper limit of the settling time: 
on a much shorter time the dust may settle down below 
the T — 1 surface of the disk, which means that settling 
can quickly affect the shape of the photospheric surface of 
the disk. In Fig. 13 the time scale is shown within which 
grains settle below the t — 1 surface. This r = 1 surface 
is computed for radially outward moving photons, using 
only the absorption opacity for silicate grains of 0.1 /im 
size at A = 0.55/im (i.e. k — 2 x lO'^cm^/g), and is shown 
as the grey curve in Figs and [2 
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Fig. 1. Left: the dust depletion height Zdopi for the disk model described in the text, for a — 0.01 and q — 1/2. This 
is the height above which virtually all grains of a certain size have been removed by settling, once the equilibrium 
settling-stirring solution is reached. The solid grey curve shows the r = 1 surface of the disk at A = 0.55/im if the dust 
has not yet settled and the opacity is dominated by grains smaller than O.l/im. Right: the time scale tsett for reaching 
such an equilibrium. This is computed from Ea. (|10|l with z = Zgctt (i-e. not with z — Zdcpi, since we are interested here 
in finding the time scale for the entire vertical equilibrium to set in). The solid grey curve is the viscous (i.e. accretion) 
time scale iaccr = /v- The dotted grey curve is the Kepler time scale tkop — ^I^k- Since the settling cannot happen 
on a time scale smaller than this ikep, we have added t^cp to the plotted time scale to ensure this minimum. The 
annotation for the black curves is the same as in the left panel. 




Fig. 2. As Fig. □ but now for a = 0.0001. 



3. 1-D vertical models: How the upper layers 
deplete 

The analytic estimates presented above give a feeling for 
the spatial and time scales on which the dust settling takes 
place. If we would only be interested in the steady state 
solution, we could also u se the analytical expressions de- 
rived by Dubrulle et al |(199ItI >. However, in order to see 
how an initially fully mixed disk evolves through settling, 
we prefer to use a numerical solution from which we can 
extract the time variations. We do this in two steps. First 
we take a close look at the local solution at a given dis- 
tance from the star. In section 01 we will then integrate 



over the entire disk in order to study the emerging SED 
and disk images. 

Using the expressions for the settling velocity and 
the turbulent diffusion coefficient presented in section |21 
it is straightforward to derive a one-dimensional verti- 
cal advection-diffusion equation for the motion of grains. 
In this one-dimensional approximation we neglect radial 
migration of th e grai ns due to gas drag or other effects 
fWeidenschilling ll977|) . Also we ignore the accretional evo- 
lution of the disk itself, or radial mixing phenomena. We 
only consider the vertical motions of the grains induced by 
settling and vertical turbulent stirring. We write the dis- 
tribution function of grains of mass m located at height 
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z above the midplane as f{m,z). It is defined such that 
/(m, z) dm dz represents the number of grains (of this size, 
at this height) per square cm of disk. 

For the gas we use again a Gaussian density profile and 
isothermal temperature profile. We assume that during 
the dust settling process this gas density and temperature 
profile do not change with time. 

The conservation equation for the dust grains is (ig- 
noring the R coordinate for compactness): 

d 

— {f{m, z)vsctt{m, 2)) = . 
oz 

This equation is valid if no dust coagulation takes place. 

We take the same stellar parameters as in Section |2 
and the same disk setup. We take our slice at i? = lOAU, 
which means that the gas surface density at that point is 
E = 13g/cm^. We simulate the settling of 0.1 /xm grains 
and adopt a — 0.0001 for the turbulence. In Fig. 0] the 
evolution of the dust settling is shown for this slice, and 
a comparison of the stationary solution (which is reached 
at about t = 3 x lO^yr) to the analytic estimates is given 
as well. 

One can see that as the grains settle, the upper parts 
of the disk are entirely depleted from dust. After only 
10'^ years the depletion above z/R = 0.30 is a factor of 
10-1° and above z/R= 0.325 it is lO'^". This means that 
long before a steady state is reached, the regions above 
z/R = 0.30 are already completely dust-free. One can see 
in Fig. 01-left that the dust that has settled forms a layer 
on top of the disk in which the dust-to-gas ratio is slightly 
above 0.01. This is simply the piling-up of the settled dust 
from higher region of the disk. 

It is interesting also to see that while the dust is still 
settling it has a very sharply defined surface. Below a cer- 
tain z the dust-to-gas ratio is about 0.01, and just a small 
step above this z the dust-to-gas ratio is almost zero. This 
sharp transition is a consequence of the fact that initially 
all grains are settling with the equilibrium settling ve- 
locity, unperturbed by stirring. Since the settling speed 
decreases towards the mid-plane, the upper-most grains 
catch up with grains just below them and therefore pro- 
duce a very sharp boundary. In the later stages, when the 
steady state has been reached, this sharp cut-off softens, 
but it should still be noticed that beyond z ~ 0.23i? the 
depletion increases very fast. 

The T = 1 surface of the disk, i.e. the disk's photo- 
sphere, for silicate grains of 0.1 /xm at t = (i.e. the 
well-mixed case without settling) lies in this example at 
about Hs/R = 0.3. Note that this is the r = 1 surface for 
the stellar radiation which enters the disk under a grazing 
angle oi (j) — 0.05. After lO'^ years the dust has settled 
well below this point, meaning that the dust will take its 
r = 1 surface along with it toward lower z (as already pre- 
dicted in Section 12). In other words: after 10^ years the 
photospheric surface height of the disk is reduced as a 



result of dust settling. The gas still remains at higher z, 
but it is the dust that carries the opacity and determines 
the disk's appearance. 

4. 1+1-D disk models: SEDs and images 

As we have seen above, dust settling can strongly affect 
the shape of the r = 1 surface Hs{R) of the disk. This 
shape determines for a large part the thermal infrared 
emission from this disk and also the appearance of im- 
ages in scattered light. If the disk's surface has a strongly 
flaring shape, one expects a spectral energy distribution 
(SED) with a strong emission, particularly at mid- and far- 
IR wavelengths. If the disk has such a shape that the outer 
regions are shadowed by the inner regions (if Hs/R de- 
creases with radius) , then the expected SED has a rather 
weak far-IR excess. The predicted images in scattered light 
for these two cases are also expected to be very different. 
It is therefore of crucial importance to investigate how the 
shape of the disk changes as a result of dust settling. 

We again start with the simple model described earlier, 
and we do a 1-D dust settling calculation at every radius, 
taking a radial grid of 20 grid points in R. In this way 
we simulate the entire disk in a 1-l-l-D manner. We take 
a = 0.00001, i.e. a factor of ten smaller than taken in 
Section 12.41 in order to accentuate the effects of settling 
on the appearance of the disk. At < =1, 10, . . . , 10^ years 
we apply a 2-D radiative transfer code (called RADMC^) 
to the given dust distribution, and compute the SED of 
the disk at these times. Also we compute the intensity 
of scattered light as a function of distance from the star, 
i.e. the brightness of an images of the disk in scattered 
light at A = 1.65/im. 

The dust temperature computed by RADMC is likely 
to differ somewhat from the dust/gas temperature T(R) 
assumed as input to the density structure of the disk. 
In principle one should iterate the computation of the 
temperature and the density structure to obtain a self- 
consistent model as we did in DD04. However, since our 
main interest is the process of dust settling rather than 
the details of disk structure, we ignore this feedback for 
the sake of simplicity. 

4.1. Model with single grain size 

The first model (model A) only has grains of O.l/im size. 
In Fig. El-left the SED as a function of time is shown. The 
SED changes only barely on time scales smaller than 10^ 
years. But on time scales longer than 10^ years the far-IR 
flux drops. Two effects play a role in reducing the far-IR 
flux. First of all, due to the setthng, the surface height of 
the disk Hg has been reduced. Since the thickness of a disk 
is directly related to the total amount of stellar radiation 
captured and reprocessed by the disk, this effect reduces 
the integrated IR emission over the entire wavelength do- 

^ www.mpa-garching.mpg.de/PUBLICATIONS/DATA/ 
radtrans / radmc / 
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0.30 0.35 



Fig. 4. The dust-to-gas ratio in a 1-D vertical slice at 10 AU in the standard disk model with a — 0.0001, assuming 
O.l/um size grains. At the start {t = 0) the dust is well- mixed with the gas, i.e. the dust-to-gas ratio is constant (we 
take it 0.01). Left: the time sequence. Right: the stationary solution compared to the analytic estimates. The vertical 
lines are the computed values of Zsctt and Zdopi (from left to right). These are defined to be solutions of Ea. (|18|l with 
^ = 1 and ^ — 100 respectively. As one can see, Zdopi is a good definition of the height above which the depletion is 
almost perfect. 



main. Chiang et al. used this fact to argue that 

dust settling may be responsible for the low far-IR flux of 
some Herbig Ae/Be star. The other, more important rea- 
son for the reduction of the far-IR flux is that due to the 
turn-over in the Hg/R beyond R = i?turn, the outer disk is 
shadowed by the inner regions of the disk. This means that 
these outer regions are irradiated only in an indirect way, 
for example by light that scattered off dust grains hovering 
slightly above the photosphere at the turn-over (i.e. self- 
shadowing) radius i?turn, or by radial radiative diffusion 
(see DD04 for a discussion of self-shadowed disks). 



Start (Only 0.1 fim grains) 




The time evolution of the scattering intensity shows 
a very interesting behavior, as is shown in Fig. [Sj-right. 
Initially, the scattering intensity changes only in the in- 
nermost regions (inward of 10 AU), from the inside out. 
After about 10^ years, the inner disk stabilizes - the steady 
state solution is reached there. By this time the outer part 
of the disk up to 100 AU starts to be affected by the set- 
tling. After lO^years, the scattering intensity is already 
reduced by more than a factor of 100, and after lO^years 
in the outermost regions, the reduction is even 10^. This 
is clear evidence of the self-shadowing effect operating in 
these regions. Without shadowing, a much more gradual 
change like the one seen in the inner regions would occur 
also here. A self-consistent computation of the disk struc- 
ture would only enhance this effect since the temperature 
in the outer disk will drop and lead to a reduced disk 
scale height. In Fig. El the self-shadowing effect is shown 
by the shape of the r = 1 surface for radially outward 
moving stellar photons. Up to 10^ years the z/R location 
of the T — 1 surface still increases for increasing R, which 
is indicative of a flaring geometry. But after 10^ years the 
T = 1 surface stays essentially at constant z/R beyond 



0.0 1 . . . . ^ 

0.1 1.0 10.0 100.0 

R [AU] 

Fig. 6. The shape of the r = 1 surface, where the optical 
depth r is at 0.55/im and is measured radially from the 
stellar surface outward (i.e. horizontally in the figure). For 
this reason the r = 1 surface is a monotonically increasing 
function of R, even when the disk fully collapses beyond 
some point. This is indeed apparent in this figure, because 
it appears if beyond 10^ and 10^ years nothing changes, 
whereas the scattering image still becomes dimmer in this 
time frame (Fig. O . 

about 7 AU, which indicates that the disk is shadowed 
beyond 7 AU. 

4.2. Model with grain size distribution 

In section 12 we have shown that the settling behavior of 
grains varies strongly with grain size. It may be possi- 
ble to represent the grain size distribution with a typi- 
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Fig. 5. The SED (left) and the 0.55/im scattering intensity (right) of the disk as a function of time for the 1+lD 
model with a single grain size of a = 0.1/im. Note that the time intervals are logarithmic (a factor of 10 apart). The 
disk is seen face-on. 



cal size for a specific application, but the correct size to 
use will depend upon the application. The intensities seen 
in scattered light is dominated by one grain size, the in- 
frared emission by another grains size, and other indica- 
tors like PAH emission will depend on yet another grain 
type. Therefore we now show the results of a calculation in 
which we include a full size distribution. The grains of each 
size settle with their own settling velocity and have their 
own Schmidt number. The bigger grains will settle much 
faster and deeper than the smaller ones, so that one ob- 
tains a different grain size distribution at different heights 
above the midplane. As our initial grain size distribution 
we ta ke an MRN distribution (Mathis, Rumpl & Nordsieck 
between grain sizes of 0.005 /im and 0.25 /im. The 
original MRN distribution is o nly well constrain ed down to 
grain sizes of about 0.02 /im l(Kim et al.l . Il99.'^ . However, 
in order to see what the influence of smaller grains would 
be we extend the size distribution down to 0.005 /im. 

The results are shown in Fig. [7| The same shadow- 
ing effects are seen as in the case of the single 0.1/im 
size grains, although there are also some differences. It 
is interesting that this shadowing in fact takes place also 
for the MRN distribution, even though the very small- 
est grains clearly stay afloat much higher up. However, 
these small grains have relatively low opacity at the stel- 
lar wavelengths: the opacity is highest for grains between 
0.1/im and l/xm. Moreover the MRN distribution goes as 
m? f(m) oc m}/^ oc a^^^, meaning that most of the mass is 
concentrated in the bigger grains. So even if the opacity of 
the grains (per gram dust) would be independent of grain 
size (up to 1 /im), still the opacity will be dominated by 
the 0.1 • • • 1/im grains. Therefore the very small grains will 
not account for much of the reprocessing of stellar radia- 
tion into the infrared. Since the albedo of the small grains 
at stellar wavelengths is virtually zero, they will not con- 
tribute at all to the scattering of starlight into the Hne of 
sight. In effect, only the grains larger than roughly 0.03/im 



will contribute significantly to the reprocessing and scat- 
tering. Therefore, even though the very small grains virtu- 
ally do not settle at all, the relevant grains (a > 0.03/im) 
do settle. 

This has two effects, both of which contribute to the 
dimming of the scattering image. First, the large grains 
which are most refiective fall below the t — 1 surface 
of the small grains (which are not reflective). The small 
grains may therefore partly shield the larger (reflective) 
grains from stellar light, thereby dimming the scattering 
image. This 'differentiation' effect is clearly seen in Fig. [3 
as an overall dimming of the disk even in regions which are 
clearly non-shadowed. This dimming is quite strong in this 
case, but for higher values of a this effect is reduced. As 
expected, such an overall dimming is absent in the single- 
grain-size case (Fig.EJ. Secondly, the t — \ surface itself is 
dragged down by settling, creating a self-shadowed geom- 
etry similar to the one for the single grain size simulation. 
The fact that even in this MRN simulation the disk be- 
comes partly self-shadowed can be seen from Fig.|Hl where 
the evolution of the t — 1 surface is shown for the MRN 
model. But it can also be seen from this figure that the 
shadowing is not as sharp as in the case of a single 0.1 
/im grain. Indeed, the weakening of the far-infrared fiux 
for the MRN case is less pronounced than in the single- 
grain model. For the scattered light intensity the effect of 
self-shadowing is still clearly seen, though the dimming by 
the 'differentiation' effect may, for these low values of a, 
be equally important for the (non-)detection of disks in 
scattered light. 

In Fig. d-left the evolution of the SED is shown for the 
MRN model. Interestingly, the mid- to far-infrared excess 
first slightly brightens (between and lO'* years). Only 
afterwards the mid- to far-infrared excess starts to decline, 
eventually becoming rather weak at 10^ • • • 10^ years. This 
initial brightening is because due to the quick settling in 
the inner region, while the outer regions take a longer time. 
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Fig. 7. As Fig. but for the disk model with an MRN grain size distribution. 
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Fig. 8. As Fig. but for the MRN model. 



the flaring angle of the surface of the disk temporarily 
increases. Only after the outer regions start to sink into 
the shadow of the intermediate disk regions the mid- to 
far-infrared excess starts to decrease significantly. 

Fig- m shows images of the disk at 0.55 /im at inclina- 
tion 45°. We truncated the disk to 100 AU so that both the 
inner regions and the outer regions can be seen clearly. 

5. Discussion 

5.1. Scaling of the results 

The analytic formulae of Section|5]and the numerical mod- 
els of Section |2| can easily be scaled to the case of fluffy 
grains. The critical parameter that goes into all these 
equations is the ratio of grain geometric cross section and 
grain mass: a/m. A compact silicate grain with a specific 
weight of X — 3.6 g/cm'^ and a size of a = 0.1/im has 
a/m = 3/(4ax) — 2 x 10^ cm^/g. A fluffy grain of size 
10/xm but the same value of a/m — 2 x 10*^ cm^/g will 
behave exactly the same as the compact 0.1 fim size grain. 
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Fig. 9. Images of the a — 0.00001 disk model with an 
MRN distribution of compact grains (see Fig.[7|), seen at 
an inclination of 45° at A = 0.55//m. For these images the 
model was truncated at 100 AU, so that both the inner 
regions (~ 25 AU) and the outer regions (w 100 AU) 
can be seen clearly. The images are at 0, 10^, 10^ and 
10^ years. The gray scale is logarithmic intensity (inverse 
video). Note that the central star is not included in these 
images (a perfect coronograph). 



5.2. The role of viscous dissipation 

Throughout the paper we have taken a very simple tem- 
perature and density structure of our disk. In reality 
the structure of an irradiated disk will be more complex 
(see e.g. Dullemond, van Zadelhoff & Natta 120021: DD04). 
Moreover, if the viscous heating by active accretion is in- 
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eluded, the midplane can become considerably hotter than 
in our simple model. While a detailed inclusion of these 
effects is beyond the scope of this paper, it is instructive 
to make some estimates. 

The midplane temperature for a disk heated 
only by viscous dissipation is given by ^^^^ = 
{27 /Q4)(aKiik/ a^mp)Y?Q,K, and the surface tempera- 
ture ("effective temperature") is given by Tv,eff = 
(8/3KiiI])i/4rv,mid, where kr is the Rosseland mean opac- 
ity of the dust-|-gas mixture (before settling). By finding 
the radius beyond which the Tv,mid and Tv^cff fall below 
the temperature given by Ea. (|2()(l . we can estimate be- 
yond which radii the viscous heating becomes unimpor- 
tant at the midplane and the surface respectively. For the 
model parameters given in Subsection l2.4l for the case of 
a = 0.01 these radii are 5 AU and 0.6 AU, for a = 0.0001 
they are 1 AU and 0.05 AU, and finally for the case of 
a = 0.00001 these are 0.5 AU and 0.02 AU respectively. 
At large radii the viscous heating is therefore not impor- 
tant, but it may be important at small radii. 

An increase in temperature due to viscous heating will 
increase the friction of the particles, and it will also in- 
crease the pressure scale height of the disk. These com- 
bined effects will increase the height of the photosphere 
after a settling equilibrium has set in. Since viscous heat- 
ing is unimportant at large radii, and therefore the in- 
crease in height of the disk takes place predominantly at 
small radii, the effect of inclusion of viscous dissipation is 
likely to increase the shadowing effect. 

5.3. Geometries 

In our earlier paper (DD04) we have demonstrated, on 
the basis of theoretical models, that disks around Herbig 
Ae/Be stars come in two main shapes: flared disks and 
self-shadowed disks. In the self-shadowed case, the inner 
rim of the disk casts a shadow over the entire disk behind 
it. The cause for the self-shadowing was a reduction of 
opacity due to grain growth. 

In the present paper we consider another process: 
dust settling. We ignore grain growth. Moreover, we fo- 
cus mainly on T Tauri stars, in which the effect of the 
puffed-up inner rim is likely to be less pronounced, and 
we ignored it here. In this setup we find geometries very 
similar, though not identical to DD04. We start out with 
a flaring disk, and after settling we may obtain a partly 
self-shadowed disk, dependent on the value of a. In con- 
trast to the self-shadowed disk solutions of DD04 our self- 
shadowed solution are only self-shadowed outside of a few 
AU, while they remain flaring inside of this radius. And 
the shadow is cast by the flaring inner part of the disk 
instead of the puffed-up inner rim. In Fig. EH the two ge- 
ometries are shown. 

For Herbig Ae/Be stars the inner rim cannot be ig- 
nored that easily. It may happen that dust settling could 
cause the entire disk to fall into the shadow of the inner 
rim, producing the fully-self-shadowed geometry presented 



in DD04, but in a different way than was shown by DD04. 
However, dust settling also affects the puffed-up inner rim, 
and this is more difficult to model in detail. The inner rim 
is much hotter than the rest of the disk. A higher abun- 
dance of free electrons would cause the magneto-rotational 
instability to operate more efhciently here and increase 
turbulent stirring. Including the inner rim into the analy- 
sis is therefore a subtle matter, which we defer to a future 
paper. 



5.4. SED type versus scattering images 

In DD04 we showed that for Herbig Ae stars, disks domi- 
nated by small grains usually have flared geometries, while 
disks in which most of the small grains have coagulated to 
big grains are usually self-shadowed. The SEDs of these 
two kinds of disks were very reminiscent of the two types 
of SEDs found among the SEDs of Herbig Ae/B e star s 
(the 'group F and 'group IF sources of Meeus et al. l200l|) . 
It was therefore suggested that group I sources are flared 
disks and group II sources are self-shadowed disks. This 
suggestion has passed several observational tests, includ- 
ing the identification of UX Orionis stars with group II 
sources (Dullemond et al.|2003), the presence or absence 
of PAH emission features (Acke et al. A&A subm.), and 
the correlation of the group I/II identification with the 
(sub-)mm slope of the source (Acke et al. subm.). The 
self-shadowed solutions in DD04 were caused by a simple 
optical depth effect: when most of the small grains are 
turned into a few big ones, the optical depth is greatly 
reduced. The disk then gradually becomes self-shadowed, 
since at large radii the optical depth is not sufficient to 
keep its surface above the shadow of the inner rim. In the 
current study we now show that also dust settling may 
turn a flared disk into a self-shadowed one. In reality, it is 
most likely a combination of settling and coagulation leads 
to the observed effects. 

It was mentioned in DD04 that self-shadowed disks 
are expected to have weak scattering images, for the ob- 
vious reasons that direct stellar radiation cannot irradiate 
the disk. But it was also shown that only for extreme pa- 
rameters the self-shadowed disk solutions that follow from 
well-mixed models are very deeply shadowed. We estimate 
that for a fully mixed model, the scattering intensity of the 
disk is reduced by typically a factor of 10. In the present 
paper we now show that dust settling can cause a much 
stronger shadowing effect. It can remove all grains capable 
of scattering stellar radiation (i.e. grains larger than about 
0.03 fim size) from the upper layers of the disk, so that if 
self-shadowing occurs in this way, the expected scattering 
images are extremely weak. Therefore, group II sources 
have weaker scattering images than group I sources. If, 
on top of that, it is found that the group II sources are 
orders of magnitude weaker, then it is clear that dust set- 
tling must also play a role. 
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Fig. 10. Pictograms of the two geometries. Left: the flaring disk we start out with. Right: after setthng the outer parts 
of the disk are self-shadowed. Note that whether this partly self-shadowed geometry is obtained depends whether the 
turbulence is low enough. 



5.5. Timescales 

The most realistic simulation in our study is the run in- 
volving an MRN grain size distribution. As shown in fig- 
ures and IHl the changes in the SED and in the scat- 
tering image brightness happen largely on a timescale of 
10^ years. Given the fact that the observed T Tauri stars 
are typically between 1 and lOMyrs old, this gives an in- 
dication that settling might be happening too fast in our 
simulation. If indeed the settling has happened already af- 
ter 10® years, one would expect most disks to show steep 
SEDs and weak scattering images. Looking at the sample 
of Herbig stars studied bv lMeeus et al.1 l)200l[l . this seems 
to be consistent in the way that the group II sources^ 
have ages larger than 1 Myr while two of the group I 
sources with available ages are only about 10^ years old. 
However, with AB Aur (about 2 Myrs) and in particular 
HD 100546 (about lOMyrs), there are two stars in the flar- 
ing disk group I which are significantly older. For T Tauri 
stars, the situation is less clear as the SEDs of T Tauri 
stars have not yet been classified in the same way as the 
Herbig stars so far. There may be observational biases in- 
volved here: The stars have been selected on IRAS colors 
which favors the detection of group I sources. If there is 
a large number of sources with strongly flaring disks at 
ages much larger than 10® years, it is possible that the 
dust settling study presented here is incomplete. A possi- 
ble solution could be larger friction between dust and gas 
in the upper disk layers, due to electrostatic interactions 
between charged grains and the ionized plasma. Another 
intriguing possibility is that even th e larger in terstellar 
grain s are already porous aggregates l)Mathis fc Whiffcn. 
Il989t IStepnik et all. l2002() — an effect which would in- 
crease their surface area and therefore their drag coupling 
to the gas. And finally, it could be possible that, for some 
unknown reason, some disks are more turbulent than oth- 
ers. The more turbulent disks then remain flaring while 



^ 51 Oph is a special case of a star with an extremely low IR 
excess and is therefore not considered here. 



the quiescent disks will become self-shadowed in a very 
short time. If this is true, then this would imply that 
group II sources are not evolutionarily linked to group 
I sources. However, since many indications point toward 
such an evolutionary link, and since there seems to be 
no good reason why some disks are more turbulent than 
others (other properties remaining the same), we consider 
this scenario less likely. 

6. Conclusion 

We find that dust settling can change the shape of the 
photosphere of protoplanetary disks on time scales of less 
than 10'^ years at small radii (< lOAU) and 10^ years 
at larger radii. On a time scale shorter than 10® years a 
equilibrium state sets in in which settling is balanced by 
turbulence. In the equilibrium state the photosphere of the 
disk may have a partly self-shadowed shape, so that the 
outer regions of the disk are not anymore directly exposed 
to the stellar radiation. They are only indirectly heated by 
infrared emission from other parts of the disk. 

Close to the star the settling in the upper layers of 
the disk proceeds much faster than at large radii, and 
an equilibrium solution is therefore reached more quickly. 
This means that the reduction in the photospheric height 
of the disk {Hs{R)) first happens at small radii, and then 
happens at progressively larger radii, until after about 10^ 
to 10® years the equilibrium state has been reached glob- 
ally. The transition from flaring disk to half-flaring-half- 
self-shadowed disk therefore happens on the latter time 
scale. 

The flaring disks (before settling) are expected to have 
a reasonably strong mid- to far-infrared excess, while the 
partly self-shadowed disks have a much weaker mid- to far- 
infrared excess. Moreover, the partly self-shadowed disks 
should become virtually undetectable in resolved scattered 
light images beyond the self-shadowing radius, whereas 
non-self-shadowed disks (early stages) should be much 
brighter. But even without self-shadowing our calculations 
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show that setthng can cause dimming of the disk in scat- 
tered hght. This is by a 'differentiation' process in which 
the larger/heavier more reflective grains settle deeper into 
the disk, leaving the smaller /lighter less reflective grains 
determining the optical appearance of the surface of the 
disk. 
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